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We present an exact calculation of the mean first-passage time to a target on the surface of a 2D or 
3D spherical domain, for a molecule alternating phases of surface diffusion on the domain boundary 
and phases of bulk diffusion. We generalize the results of [1 and consider a biased diffusion in 
a general annulus with an arbitrary number of regularly spaced targets on a partially reflecting 
surface. The presented approach is based on an integral equation which can be solved analytically. 
Numerically validated approximation schemes, which provide more tractable expressions of the mean 
first-passage time are also proposed. In the framework of this minimal model of surface- mediated 
reactions, we show analytically that the mean reaction time can be minimized as a function of the 
desorption rate from the surface. 

PACS numbers: 



I. INTRODUCTION 



Reaction kinetics in confined systems where a small number of react ants are involved, such as porous catalysts and 
living cells, can be limited by the time needed for molecules to meet and react [2j|3]. This time is known in random 
walk theory as a first-passage time (FPT) [4-7 . For the specific case of biochemical reactions in living cells, these 
general considerations have to incorporate two important features. First, while passive diffusion is the dominant mode 
of transport in chemical systems, active transport has been shown to play a prominent role in living cells [8]. As 
a matter of fact, various motor proteins such as kinesin and myosin are able to convert the chemical fuel provided 
by ATP into mechanical work by interacting with the filaments of the cytoskeleton. Many macromolecules or larger 
cellular organelles such as vesicles, lysosomes or mitochondria, can randomly bind and unbind to these motors [9HTT]. 
As a result, the overall transport of such tracers in the cell can be described in a first approximation as diffusion in 
a force field [12] . Second, reactions in confined domains like cells generally involve surface- mediated diffusion that 
combines bulk transport and surface diffusion due to non-specific interactions with the domain boundary (e.g. cell 
membrane) [13-18]. Such two-state paths and the corresponding first-passage properties have been studied in the 
broader context of intermittent search strategies [19-22 under the hypothesis that the times spent in each state 
(surface and bulk) are controlled by an internal clock independent of any geometrical parameter. In most cases, the 
sojourn times in each state have been assumed to be exponentially distributed [21 , with the notable exception of 
Levy [23] and deterministic laws [24, 25] . However, in the case of interfacial reactions, for which molecules react on 
target sites located on the surface of the confining domain, the time spent in a bulk excursion is controlled by the 
statistics of return to the surface and therefore by the geometry of the confining domain [26H30] . Hence this return 
time is not an external parameter but is generated by the very dynamics of the diffusing molecule in confinement. 

Recently, such coupling of the intermittent dynamics to the geometry of the confinement has been explicitly taken 
into account |T| I3TH35] . Exact calculations of the mean FPT to a target on the surface of a 2D or 3D spherical domain, 
for a molecule performing surface-mediated diffusion, have been proposed [TJ[33]. However, these works have been 
limited to passive transport alone. The present article develops a general theoretical framework which in particular 
allows one to incorporate the effect of active transport on surface-mediated diffusion. 

More precisely, we extend the results of [1, 33 in four directions: (i) we consider the general case of an imperfect 
adsorption step, so that the molecule can bounce several times before being adsorbed to the confining surface [36ti44] ; 
(ii) the geometry adopted is a general annulus, whose either interior, or exterior boundary is purely reflecting; (iii) 
we take into account the effect of an exterior radial force field which, for instance, can schematically mimic the effect 
of active transport; (iv) we consider the possibility of having an arbitrary number of regularly spaced targets on the 
surface. Relying on an integral equation approach, we provide an exact solution for the mean FPT, both for 2D an 3D 
spherical domains, and for any spherical target size. We also develop approximation schemes, numerically validated, 
that provide more tractable expressions of the mean first passage time (MFPT). 

The article is organized as follows. In Sec. [TTJ we define the model under study; in Sec. |III[ we show that the MFPT 
can be determined by solving coupled partial differential equations that can actually be converted into a single integral 
equation. We then provide an exact solution of this integral equation, as well as an approximate, more tractable, 
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solution. In Sec. IV, we give fully explicit expressions of the MFPT by applying this general formalism to particular 
cases, representative of the four aforementioned extensions. 



II. THE MODEL 

The surface-mediated process under study is illustrated in Fig. [I] We consider a molecule diffusing in the volume 
S between two concentric spheres of radii R and R c . The molecule alternates phases of bulk diffusion (with diffusion 
coefficient D2) and phases of surface diffusion on the boundary of the sphere of radius R (with diffusion coefficient 
Di) which contains a target. The target is defined in 2D by the arc G [— e, e], and in 3D by the region of the sphere 
such that 6 G [0, e] where 6 is in this case the elevation angle in spherical coordinates. Note that as soon as e ^ 0, the 
target can be reached both by surface and bulk diffusion. 

In the following, the case R > R c will be called an exit problem and the case R c > R an entrance problem (Fig. [T]). 
In 3D, the entrance problem can account for the time needed for a virus initially in the cell (the sphere of radius R c ) 
to get into the nucleus (the sphere of radius R) through a single nuclear pore (the target) in the presence of diffusion 
on the nuclear membrane. In turn, the exit problem in 3D may describe macromolecules searching an exit from the 
cell through a channel (or channels) in the cellular membrane. In that case, the surface of the nucleus is considered 
as purely reflecting. The 2D geometry could correspond to cells that are confined, as realized in vitro for example in 

In this model, a molecule hitting the sphere of radius R c is immediately reflected. In contrast, when the molecule 
reaches the sphere of radius R, which contains the target, it is imperfectly adsorbed: the molecule hitting the boundary 
at r = (R,6), 6 G [0,7r] is at random either adsorbed to the sphere of radius R or reflected back in the bulk. The 
quantity k which describes the rate of adsorption is more precisely defined through the radiative boundary condition 
Eq. Q (see also Eq. ( |A4[ ) of the discrete lattice approach discussed in Appendix|A|. In particular, k = 00 corresponds 
to a perfectly adsorbing boundary and k = to a perfectly reflecting boundary. Notice that for finite /c, molecules 
hitting the target from the bulk can be reflected. 

The time spent during each surface exploration on the sphere of radius R is assumed to follow an exponential 
law with desorption rate A. At each desorption event, the molecule is assumed to be ejected from the surface point 
r = (R, 0) to the bulk point r = [R — a, 6). In what follows, a can be positive or negative: a > for the exit problem 
(R > i? c ), and a < for the entrance problem (R < R c ). Although formulated for any value of the parameter a such 
that \a\ < \R — R c \ (to ensure that the particle remains inside the domain after reflection), in most physical situations 
of interest \a\ is much smaller than R. Note finally that a non zero ejection distance a is required in the limit of 
perfect adsorption k = 00, otherwise the diffusing molecule would be instantaneously re- adsorbed on the surface. 
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FIG. 1: Model - Left: Static picture of the entrance problem in 3D - Right: Dynamic picture of the exit problem in 2D. The 
green sphere stands for the diffusing molecule and the red sector stands for the target. 
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III. GENERAL SOLUTION 
A. Basic equations 

For the process defined above, the mean first-passage time (MFPT) satisfies the following backward equations 

T ^^ e t 1 {e) + \{t 2 {R-a,0)-h{9)) = -1 (e<0<w), (1) 

D 2 (A r + ^ d r + ^f) t 2 (r,e) = -1 ((r,e)eS), (2) 

where: (i) t\(0) stands for the MFPT for a molecule initially on the sphere of radius R at angle and ^(r, 0) stands 
for the MFPT for a molecule initially at a bulk point (r, 0) within the annulus S = (R C ,R) x [0, 7r]; note that, due 
to the symmetry ti{6) = ti(—0), in 2D 6 can be restricted to [0,7r]; (ii) the radial and angular Laplace operators are 
respectively 



dr 2 



■ — , A^ = (sin(9) 2 - d ^(sin(9) G 



and d stands for the space dimension (in practice, d will be taken equal to 2 or 3); (hi) v(r) is the radial velocity of 
the molecule resulting from an external force. 

In Eqs. ([lj[2|, the first terms of the left hand side account for diffusion respectively on the surface and in the bulk, 
while the second term of Eq. describes desorption events. These equations have to be completed by boundary 
conditions: 

(i) reflecting boundary condition on the sphere of radius R c > 

^ =0 (0<0<tt) (3) 

or \ T =(R C ,6) 

(note that this condition holds even in the presence of the velocity field v(r), see e.g. [45]): 

(ii) radiative boundary condition 

^ =k{t 1 (0)-t 2 (R,0)} (0<#<^), (4) 

dr \r=(R,e) 

which describes the partial adsorption events on the sphere of radius R (see Appendix [A] for justification of this 
boundary condition). For the exit problem (R > R c ), the radial axis points towards the surface and k > 0, while for 
the entrance problem (R < R c ), the radial axis points outwards the surface and k < 0. Finally, the limit k = ±oo 
describes the perfect adsorption for which the above condition reads as t\{6) = t2(R,0). 

(iii) Dirichlet boundary condition 

*i(0) = O (0<6><e), (5) 

which expresses that the target is an absorbing zone (the search process is stopped on the target). 
In what follows we will use two dimensionless quantities 

x = l- a/R, (6) 

UJ = R^/D[, (7) 

and the operator L acting on a function / as 

(L/)(r) = /(r-a)-/(r)-i&/(r). (8) 

B. General integral equation 

We generalize the approach presented in [1 and show that the coupled Eqs. ([ll [2]) with the boundary conditions 



3l|5]) lead to the integral equation (22) for t\ only. 
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The starting point is a Fourier decomposition of ^. Eq. Q is easily shown to be satisfied by 

^ oo oo 

t 2 (r, 0)=a + —f(r) + A) /o(r) + ^ a n f n (r)V n (0) + ^ «-n/-n(r)^ n (0), (9) 



with coefficients a n to be determined, and 

(i) f(r) is a rotation-invariant solution of Eq. Q regular at r = 0, i.e. 



A r + ^ a r ) />) = -i, (io) 

the choice of f(r) being up to an additive constant; 

(ii) fo(r) is a non-constant solution of the homogeneous equation 

r + v ~ih dr ) /o(r) = °' 

the choice of fo(r) being up to an additive constant and a multiplicative prefactor. It can be shown that fo(r) 
necessarily diverges at r = in our cases of interest; 

(hi) the set of functions {/n(0? Ki(#)}nez is an eigenbasis of the homogeneous equation associated to Eq. (J2|: 

- A e V n (0) = p n V n {6) (n > 0), (12) 

'•^A r +^9r)/ n W = p H /nW (neZ), (13) 

with V- n {9) = V n {9) due to the reflection symmetry, and 

(14 » 

n(n + l) (a = 3). 



r 



We set 



1 (n = 0) (d = 2) 

V n (9) = { [V2cos(n6) (n > 0) 1 ; ' (15) 

y2n + 1 P n (cos0) (n>0) (d = 3), 

where P n (z) are Legendre polynomials. In tur n, t he functions / n (0 which depend on the velocity field u(r), will be 
determined individually case by case (see Sec. IV). 
In the following, we will use two inner products: 



(/,<?) -»• (f\g)= I' f(9)g(0W d (9), 

JO 

(f,9) -»• (f\9)e = f f(9)g(9W d (9), 



where d/j,d(0) are the measures in polar (d — 2) and spherical coordinates (d = 3): 

d9 sin /9 

= -, d^(9) = —d9. (16) 

7T Z 

With these definitions, the eigenvectors V n {0) are orthonormal 

(K|V m ) = 5 nm . (17) 

We now use the boundary conditions ([3][5| to determine the coefficients {a n } n defining ^(r, 0) in Eq. 
(i) The reflecting boundary condition (|3| reads 

A) d rfo(r)\R c + -p-drf(r)\R c + ( a n d rfn + ^-nd r f-n)\R c V n (0) = 0, (18) 

2 n=l 
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which, using the orthogonality in Eq. (17), leads to the following relations 

1 



_ d r f(r) 

D 2 Ur/oW 



a n d r f n (r)\ r=Rc = -a- n d r f - n {r)\ r=Rc . 



(19) 



\r=R c 



Note that, in the case R c = 0, the first condition reads /3q = 0. Indeed, if (3q were non zero, the MFPT of a molecule 
initially at the origin would diverge (by definition of the function /o). 

(ii) Substituting Eq. Q into the radiative boundary condition Eq. Q, projecting it onto the basis V n (9) and using 
Eq. (19), we obtain two supplementary conditions: 



a 



1 I d r f(r) 
D 2 Ur/o(r) 



\r=R c 



f n (R) + T d r f n (R) 



MR) + ^drfo(R) 



9 r fn(r) 



1 

D~2 



f(R) + -d r f(R) ) = 



d r f-„(r)J lr 



f- n (R) + ^d r f- n {R) 



= (ti\V n ) (n>0). 



(20) 



(21) 



On the other hand, the radiative boundary condition in Eq. Q and the operator L defined in Eq. ^ allow one to 
rewrite Eq. ([!]) as 



-AehiO) = ^-+u 2 (Lt 2 )(R) (e<0< tt), 



which becomes, using Eqs. ([9j[l9j[20j[2l}, 

oo 

-AehiO) = u 2 T + u 2 J2 x n(ti\V n )V n (9) (e < < tt), 



71=1 



where 



Vd 



1 v± 

( d r f(r) ^ 
Ur/o(r) 



Lfo(R) + Lf(R), 



\r=R c 



fn(R) + \drUR) - (ifiSl)) (f-n(R) + R/-n(i?)) 



(n > 1). 



(22) 

(23) 
(24) 

(25) 



In Appendix |B] we identify the quantity r] ( i/D 2 as the mean first passage time on the sphere of radius R for a molecule 
initially at r = R — a. Thus the time T is the sum of a mean exploration time r]d/ 1 D 2 and a mean "exploitation" time 

1/A - 

(iii) The absorbing boundary condition (I5| and the relation t'^ir) = which comes from the invariance of t\ under 
the symmetry 9 — >• 2ir — (9, lead after integration of Eq. (22) to 



h(0) 



WT g € (9) + co 2 J2Z i f^<K|ri)e{%(0) - V n (e)} (e < 9 < tt), 
10 (O<0<e), 



where p n is defined in Eq. ( 14 ) and g e (9) is the solution of the problem: 

A e g e (6) = -1, with g e (e) = and g' € (tt) = 0. 



(26) 



(27) 



Note that R 2 g e (9)/Di represents the MFPT to the target when A = 0, i.e. in absence of desorption events, hence 
g e (9) is well known: 



9e(0) = { 



-(9-e)(27r-e-9) (d = 2), 



In 



1 - cos(<9) 
1 — cos(e) 



(d = 3). 



(28) 



Equivalently, Eq.(26) reads 

m = < 



X 

g e (0) + to 2 V ^(V n \^) € {V n {6) - V n (e)} (e < < tt), 
(0 < < e), 



where tp(9) = ti{9) / {uj 2 T) is a dimensionless function. 



(29) 



C. Exact solution 

The function i/j(0) can be developed on the basis of functions {V n {6) — V n (e)} n , 

oo 

V>(0) = g € (6) + J2 d n {V n (6) - V n (e)} (e < < n), 



with coefficients {d n }n>i to be determined. Due to Eq. (26), the vector d = {<i n } n >i is a solution of the equation 

oo oo / oo \ 

dn{V n (9) - V n (e)} = u 2 Y / [ U n+Yl Qn, m d m {V n (0) - V n (e)}, (30) 



n=l n=l \ m=l 

where we have defined the vectors U and £ by their n-th coordinates: 



Pn 



£», dn = Pn(9e(0)\V n (6)) e (n > 1), 



and the matrices Q and I e by their elements: 

X, 



Q„, m = h{n,m), I e (n,m) = (V n (0)\V m (0) - V m (e)) t (ro > 1, n > 1). 

Pn 



(31) 



(32) 



As Eq. (30) is satisfied for all G (e, 7r), the coefficients d n can be found by inverting the underlying matrix equation 



as 



d n 



(I-uj 2 Q) 1 U 



The MFPT t\{9) can be explicitly rewritten as 

oo 

g e (e) + J2d n {V n (6)-V n (e)} 



ti(0) 



tu 2 T 




n=l 



(e < < tt), 
(0 < < e). 



(33) 



(34) 



The averaged MFPT {t\) which is defined by averaging over a uniform distribution of the starting point, is then easily 
obtained as 



PIT / OO \ 

<tl> = ^ h{0)dfX d {0) = L0 2 T ( (g e \l) e + J2^nZn\ , 



(35) 



where we have used the following relation 

(V n (0) - V n (e)\l) e = -(V n (0) - V n (e)\A e9e (0)) e = p n (V n (0)\g e (0)) e = £„ 



Finally, the MFPT t2(r,6) is given by Eq. (|9|, in which the coefficients fio and a n are obtained from Eqs. (19, 20 
2li 



r=R c 



fo(r)-fo(R-a) 



+ 



OO /- 

ra=l ^ 



(r) 



&rfn 
&rf—n J r 



f-n(r) , 



(36) 
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Expressions 



2D 



3D 



V n (0) 



1 (n = 0) 

\/2cos(n6>) (n > 0) 



V2n+ 1 P n (cos ^ 



/On 



i(n+ 1) 



smOdO/2 



7 l-cos(6>)~ Y~ 
^ l-cos(e) J 



§(0-e)(27r-e-0) 



In 



0/e|l>e 



log i 



(<?e|K)e 



^2 



^2-{(7r — e) cos(ne) + sin(ne)/n} 



2n + l / p / \ uP n (u)-P n _ 1 (u) F n (u) + 1 
2 ^ n+1 t 2n+l 



7 e (n, n), n > 1 



1 , i sin 2ne \ 

— 7T — e H ~ — 

7r V 2n / 



7 x sinfme) 7 x sin(ne) 
2 cos(ne) - --cos(me) - - 2 

— 1IL —9. — —9. — TU 



I e (n, m / n) 
m, n > 1 



V2n+lV2m+l (n-m)uP m (n)P TL (tt) + (m+l)P m (u)P n _ 1 (u) - (n + l)P n (u)P m _ 1 Q) 
2 771 



(n+l)[m(m + l)-n(n+l)] 

where ^ = cose and function F n (u) is given in Appendix |c| 



TABLE I: Summary of formulas for computing the vector £ and the matrix Q in Eqs. (31 32) that determine the coefficients 
d n according to Eq. (33). 



with 



Tp n d n 



L ^( R )-(wh) r=R Lf-n(R) 



(n > 1). 



Table [j] summarizes the quantities which are involved in Eqs. p4j [35j [36) and independent of the detail of the 
radial bulk dynamics. In turn, the quantities T, rjd and X n are expressed by Eqs. p3{ [24| [25] ) through the functions 



/, fo and / n and thus depend on the specific dynamics in the bulk phase and will be discussed in Sec. 
particular examples. 



IV 



for several 



A numerical implementation of the exact solutions in Eqs. (34 1 35 , 36) requires a truncation of the infinite- 
dimensional matrix Q to a finite size N x N. After a direct numerical inversion of the truncated matrix (J — uj 2 Q) 
in Eq. (33), the MFPTs from Eqs. p4{ [35] [36]) are approximated by truncated series (with N terms). We checked 
numerically that the truncation errors decay very rapidly with TV. In a typical case of moderate u < 100, the results 
with N = 100 and N = 200 are barely distinguishable. In turn, larger values of uj (or A) may require larger truncation 
sizes. In the following examples, we used N = 200. In spite of the truncation, we will refer to the results obtained by 
this numerical procedure as exact solutions, as their accuracy can be arbitrarily improved by increasing the truncation 
size N. These exact solutions will be confronted to approximate and perturbative solutions described in the next 
subsections. 



D. Are bulk excursions beneficial? 



Before considering these perturbative and approximate solutions, we address the important issue of determining 
whether bulk excursions are beneficial for the search. This question can be answered by studying the sign of the 
derivative of (ti) with respect to A at A = 0. In terms of Q = —QR 2 /Di, the MFPT from Eq. (35) reads 



D 2 



(t 1 ) = -z(l + \ri d /D 2 ) 



Di 
R 2 



<<7e|l)e 



£-A(/ + AQ) -1 f7 



The derivative of (t\) with respect to A is 



If the derivative is negative at A = 0, i.e. 



D 2 R? 



' (r ? - 1 +2A)J + A 2 Q ^ 
(J + XQ) 2 . 



Pi 

D 2 



< 



R 2 (±U) 

(&|l)e Vd 



(37) 



(38) 



(39) 
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bulk excursions are beneficial for the search. Explicitly, the critical ratio of the bulk-to-surface diffusion coefficients, 
below which bulk excursions are beneficial, is 



Pic 



%(ge|l) 

R 2 



oo 



X n (g € \V n )l 



(40) 



E. Perturbative solution (small e expansion) 



While Eq. (34) for t\ is exact, it is not fully explicit since it requires either the inversion of the (infinite-dimensional) 
matrix I — uj 2 Q, or the calculation of all the powers of Q. In this section, we give the first terms of a small e expansion 
of the MFPT, while in the next one we provide an approximate solution that improves in practice the range of validity 
of this perturbative solution. Both solutions rely on the orthogonality of functions V n in the small target size limit 
e —> 0, which implies that the matrix Q is diagonal in this limit. 

More precisely, as V n (0) = V n (—9), necessarily doV n (0) = 0, so that for e close to zero and for all G [0, e], one has: 
V n (0) = V^(0) + 0{9 2 ). As a consequence, the function J e (n, m) introduced in Eq. (32), reads for all m, n > 1 (see 
also Appendix |C|) 

J e (n, m) = (V n (0)\Vm(0) - V m (e)) e = (V n (0)\V m (0)) - V m (e)(V n (0)\l) + 0(e 3 ) = S nm + 0(e 3 ). (41) 
The first terms of a small e expansion of the MFPT can then be exactly calculated. Relying on the expansion Eq. 



(41), one can replace I e (n, n) by 1 to get in 2D 

^ = ( 4 + 2 ^ E 2/ 2 Xn 2v^ - ™ + ( 1 - 2 ^ E 2 ^2^ l ^ + °^)' ( 42 ) 

cj 2 T \ 3 ^ n 2 (n 2 - uj 2 X n ) \ ^ n 2 - u 2 X n v ' v 7 

\ n=l v 7 / \ n=l / 

and in 3D 

^ = -2 ln(e/2) - f 1 + c^ 2 V ^ ^t 1 ^ + 0(e 2 ). (43) 

cu 2 T V ^n(n + l)(o; 2 X n -n(n + l))y V ; V } 

The comparison of the perturbative solutions to the exact and approximate ones is presented in Figs. [2j [3J [5| [7] and 
it is discussed below. 



F. Approximate solution 



As mentioned above, we now provide an approximate solution that improves in practice the range of validity of the 
perturbative solution. This approximation relies on the fact that, due to Eq. (41), the matrix Q defined in Eq. (32) 
reads 



Qnm — 3mnQnn H~ 0{t )• 

Keeping only the leading term of this expansion, one gets 

uJ 2 U n 



(n > 1). 



From Eqs. ( |3Tj [35j [45] ) we then obtain the following approximation for the search time: 

X n (9e\V n ) 2 e 



<*1 



j 2 T 



(9e\l) 



j 2 ^I e (n,n) 



(44) 



(45) 



(46) 



Note that this expression is fully explicit as soon as the functions /, fo and f n defined in Eqs ( [Iq| [TT| [l3| are 
determined. In Section [IV| we will consider particular examples and write these functions explicitly. As we will show 
numerically, this approximation of ti, which was derived for small e, is in an excellent quantitative agreement with 
the exact expression for a wide range of parameters and even for large targets (see Figs. [2j [3j [5| [7]) . 
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Case 


Quantity 


2D 


3D 


No bias (V = 0) 


/ 


-r 2 /4 




-r 2 /6 




fo 


lnr 




i?/r 




fn 


r n 




r n 




f-n 






r ~n-l 


Velocity field: 


f 


-r 2 /(2(2- M )) 




-r 2 /(2(3- M )) 


v(r) f 


fo 


[{r/RY - l]//i 




(rWV(l-/x) 




fn 


r /x/2+ 7ri 




r (/x-l)/2+7n 




f-n 


r /x/2- 7r , 




r ( M -l)/2- 7n 






(7n = V^ 2 +M 2 /4) 


{in 


= V^^+^ + ^-l) 2 /^ 


Sector of angle 


f 


-r 2 /4 




-r 2 /6 


(no bias, V = 0) 


fo 


lnr 




i?/r 




fn 






r -(l/2)+ 7n 




f-n 






r -(l/2)- 7n 








{in 


= Vn(n + l)(7r/0) 2 + l/4) 



TABLE II: Functions /, fo and f n for several particular cases in 2D and 3D (see also Appendix |d| 



IV. PARTICULAR CASES 



We now show how the above theoretical approach can be applied to various important examples. The only quantities 



oacn can be applied to 7 

needed to obtain fully explicit expressions of Eqs. (|35j |46j |42j |43| are the functions /, fo and f n defined in Eqs. 
(10 1 11 , 13) which are involved in the definitions of the quantities T and X n according to Eqs. (23 1 25). These 



quantities are listed in Table [TT] for the representative cases discussed in this section. Throughout in this section, all 
the quantities {R, R C1 a, e, A, Di, D 2l fc, (ti)) are written in dimensionless units. The physical units can be easily 
retrieved from the definitions of these quantities. 



A. Zero bias (V = 0) 

1. Exit problem for a perfect adsorption 



In the case of the exit problem with R c = 0, perfect adsorption (k = oo) and no bias, the formula (46) reproduces 
the results of [1 . The coefficient r]d/D2 is the mean first passage time to the sphere of radius R, starting from 
r = R — a, 



rjd cl{2R - a) 
D 2 ~ 2d 



From the expressions for the quantities T and X n , 

~ i R 2 



X 2dD 2 

we retrieve the approximate expressions for the MFPT in 2D 

co 2 T 



(1 - x 2 ), X n = x n - 1, 



(h) 



1 3 2uj 2 ^ x n — 1 {{n — e) cos(ne) + sin(ne)/n) 



TT 



71=1 



1 



sin 2ri 



and in 3D: 



In 



1 — cos e 



1 + cos e 



c^y> ( x n -l)(2n + l) 



((l + ^?f)Pn(cOS6) 



i (cos e) 



n+1 



n=l 



n 2 (n + l) 2 



1 a£ (s"-l)(2n+l) f ( x 
1 2 n(n+l) 1 e\ n i n ) 



(47) 
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We emphasize that bulk excursions can be beneficial for the MFPT even for the bulk diffusion coefficient D2 smaller 
than the surface diffusion coefficient D\ [1 . This can be understood qualitatively by the fact that bulk diffusion induces 
flights towards remote and unvisited regions of the sphere r = R. These long-range hops can diminish the time for 
target encounter (provided that the time spent in the bulk phase is not too large). 



2. Exit time for a partial adsorption 



We now give an explicit expression of the results (|46|) and (40) for a 2D exit problem with R c = and with an 



imperfect adsorption on the sphere of radius R. Using the expressions from Table [01 the coefficients r]d and X n are 



Vd 



2d 



2 

kR 



1 



kR 



kR 



Thus the approximate MFPT in 2D reads 



<*i> 



j 2 T 



2oj z 

7T 



E 



kR 



((tt — e) cos(ne) + sin(ne)/n) 2 



(1 + 0)1 



* « 2 (i+ra) 



sin 2ne ^ 



(48) 



(49) 



and the critical ratio of the bulk-to-surface diffusion coefficients in Eq. (|40j), below which bulk excursions are beneficial, 
takes the form 



2 \ 7r(7r - e) 3 
kR) 24 



E 



kR ((tt - e) cos(ne) + sin(ne)/n) 2 



(1 + 0) 



.n=i ' 6 V" 1 fci*; 



(50) 



Similarly, one can write explicit formulas in 3D. 

The MFPT as a function of the desorption rate A is shown on Fig. [2] for different values of the bulk diffusion 
coefficient D2 and the target sizes e = 0.01 and e = 0.1, both in two and three dimensions. One can see that the 
approximate solution ( [49] ) (shown by circles) accurately follows the exact solution (shown by lines) for a wide range 
of parameters. In turn, the perturbative solutions in Eqs. (42, 43| (shown by pluses) are accurate for small e = 0.01 
but they deviate from the exact ones for larger e = 0.1. 

The quality of the approximate and perturbative solutions can also be analyzed on Fig. [3] which shows the MFPT 
as a function of the target size e (with a moderate value A = 10). Once again, the approximate solution is very 
accurate for the whole range of e, with a notable deviation only at e close to 1. The perturbative solution starts to 
deviate for e > 0.1 (as the desorption rate A appears in the coefficients of the perturbative series, the validity range 
would of course depend on A used). 

The situation of quasi-perfect adsorption (kR ^> 1) can be shown to be asymptotically equivalent with the case of 
short ejection distance (a/R <C 1), as illustrated on Fig. [4j 



3. Reflecting boundary and entrance time 



Now we provide an explicit form for Eqs. (46, 401) in the presence of a perfectly reflecting sphere of radius R c . We 



recall that the case R c > R (resp. R c < R) is called an entrance (resp. exit) problem. 
Using the expressions from Table \H\ the coefficients r]d and X n can be written as 



I? 
4 

6 



2 

kR 
2 

kR 



R 2 C 
2 

3R 



ln(x) 



1 

kR 
1 

" kR 



(51) 
(52) 



and 
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FIG. 2: MFPT (ti) as a function of the desorption rate A for domains with partial adsorption k = 1: comparison between 
the exact solution (lines), approximate solution (circles) and perturbative solution (pluses) for 2D (left) and 3D (right), with 
e = 0.01 (top) and e = 0.1 (bottom). The other parameters are: R = 1, D\ = 1, a = 0.01, no bias (V = 0), and D 2 takes three 
values 0.5, 1 and 5 (the truncation size is N = 200). 




FIG. 3: MFPT (ti) as a function of the target size e for domains with partial adsorption k — 1: Comparison between the 
exact solution (lines), approximate solution (circles) and perturbative solution (pluses) for 2D (left) and 3D (right). The other 
parameters are: R = 1, D± = 1, a = 0.01, A = 10, no bias (V = 0), and D2 takes three values 0.5, 1 and 5 (the truncation size 
is N = 200). 



It is worth noting an interesting dependence of (ti) on the radius R c when R c and a are both small. One finds in 
2D 

8(h) d'ih) _(D 2 7r 2 \ AXaR 



dR c \ Rc=0 ' dRl lRc=0 \Di 24J DM 
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(35) 



FIG. 4: MFPT (ti) computed through Eq 
parameters k and a, with R c = (left) and i? ( 
for small a. If both values of a and 1/fc are close to zero, the MFPT (green line) tends to be a constant which is equal to the 
MFPT on a segment of length 2tt. Here d = 2, R = 1, Di = 1, D 2 = 5, e = 0, no bias (V = 0) (the truncation size is N = 200). 



as a function of the desorption rate A for several combinations of the 
yz > R = 1 (right). The relation between 1/k and a is asymptotically valid 



and in 3D, 

dR c{Rc=0 dRl lRc=0 9i?3 |fic=o 27 L -"' V "'- ^ 1 J8D 1 D 2 

In 2D, as long as D2/D1 > 7r 2 /24 « 0.411 introducing a reflecting sphere of small radius i? c increases the search time. 
This can be understood as follow: increasing R c <C R increases the duration of flights between remote and unvisited 
regions of the sphere r = R, as these flights have to circumvent an obstacle at r = R c . These long-range flights can 
reduce the search time only if they are not too time costly, hence the condition on D2 > £>2c- The critical diffusion 
coefficient D2 C increases with R C <R (Fig. |6|. 



B. Case of a 1/r velocity field 



We now examine the case of a radial 1/r velocity field v(r) characterized by a dimensionless parameter fi\ 



v { r ) = Y~ r - 



Substituting the functions /, Jo and f n from Table [ll| into Eq. (23), we can write the coefficients rjd and X n as 



where 



m 



R 2 



2(2 - fj) 

R 2 
2(3 - /i) 



l-x z 



1-x 2 



2 \ 2 (R x 



kRJ n\R 



kR 



fi-1 \ R 



2-/i 



R r 



—) 

kRJ 



3—fi 



X n — 



7n +70 + 7n +70 f Rc 

kR 7 n - 70 V R 



2 7 n 



kR 

In ~ 70 



7« + 7o , 7» +7o / 
fci? 7n - 70 \ -R 



2 7n 



7n -70 



7n 



Vn 2 + ^ 2 /4 (d = 2) 

v/n(n + 1) + (/x - l) 2 /4 (d = 3) 



kR 



(n>l), 



(57) 

(58) 
(59) 

(60) 



(61) 



and 70 = /i/2 in 2D and 70 = (/i — l)/2 in 3D. Note that in the limit /i = 0, one gets 7 n = n (n > 0) in 2D, and 

70 = —1/2 and j n = n + 1/2 in 3D, so that the above results are reduced to the previous case. The case \i = d has 

2 

to be considered separately because f(r) = ^-(1 — 2 In r) in both 2D and 3D. 
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FIG. 5: MFPT (ti) as a function of the desorption rate A for an annulus with the inner radius R c = 0.5 and the outer radius 
R = 1: comparison between the exact solution (lines), approximate solution (circles) and perturbative solution (pluses) for 2D 
(left) and 3D (right), with e = 0.01 (top) and e = 0.1 (bottom). The other parameters are: D\ — 1, a = 0.01, k = oo, no bias 
(V = 0), and D2 takes three values 0.5, 1 and 5 (the truncation size is N = 200). 
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FIG. 6: Critical ratio of the b ulk- to-surface diffusion coefficients Dzc/Di in 2D, with perfect adsorption and no bias (k — 00, 
V = 0) computed through Eq. (40) as a function of the target size e for different values of a and R c : the exit problem (R c < 1) 
on the left and the entrance problem (R c > 1) on the right (the truncation size is N = 200). 



The same expression for rjd stands in the cases \i = in 2D and \i = 1 in 3D: 

When R c = and \i > d, the MFPT to the sphere 77^/^2 diverges, which causes the critical bulk diffusion coefficient 
D 2c to diverge. 

Figure [7] shows the MFPT (h) as a function of the desorption rate A in the presence of a 1/r velocity field. As 
earlier, the exact, approximate and perturbative solutions are in an excellent agreement for a wide range of parameters. 
Figure [8] shows a similar dependence for different field intensities \i (if \i > 0, the velocity field points towards the 
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FIG. 7: MFPT (ti) as a function of the desorption rate A in the presence of a 1/r velocity field: comparison between the exact 
solution (lines), approximate solution (circles) and perturbative solution (pluses) for 2D with \i — 1 (left) and 3D with /i = 2 
(right), with e = 0.01 (top) and e = 0.1 (bottom). The other parameters are: D\ = 1, a = 0.01, k = oo, and D2 takes three 
values 0.5, 1 and 5 (the truncation size is N = 200). 



origin, while /i < means that the velocity field points towards the exterior). For R c < R (resp. R c > R), for a fixed 
A the search is on average faster as \i is more negative (resp. positive). Finally, in Fig. [91 the critical ratio of the 
bulk-to-surface diffusion coefficients is shown as a function of the target size, both in two and three dimensions. The 
dependence on the field intensity \i is stronger in 2D than in 3D. 

For R < R C1 large absolute values of the drift coefficient increase (£1) and D2 C as (i) a strong outward drift 
(|/i| ^> |/i c |) diminishes the probability for fast relocation through the central region; (ii) a strong inward drift 
(/i ^> |/i c |) traps the diffusing molecule in the central region and increases 77^/^2, the MFPT to the surface r = R 
after ejection. 

Although we derived the formulas for both 2D and 3D cases, the 1/r velocity field is mainly relevant in two 
dimensions as being a potential field. In three dimensions, the potential field exhibits ljr 2 dependence. This case, 
as well as many others, can be treated by our theoretical approach after solving Eqs. (10 1 11 |l3) for the functions 
fo( r ), f( r ) an d f n ( r )- This is a classical problem in mathematical physics. For instance, the aforementioned velocity 
field 1/r 2 in three dimensions involves hypergeometric functions, as shown in Appendix |d| 



C. Circular and spherical sectors 



The above approach can also be applied for investigating the MFPTs in circular and spherical sectors of a given 
angle <p (Fig. 10). In most biological situation such as viral trafficking, <j> <C tt but the arguments presented here 
stand for arbitrary <j). For this purpose, the angular basis functions V n {6) can be rescaled by the factor </>/7r: 



1 



V n (9) 



y/2 cos(n#7r/(/>) 
y2n + 1 P n (cos (07r/0)) 



(n = 0), 
(n > 0) 
(n > 0) 



2) , 

3) , 



(62) 



15 



3.5 
3 
2.5 

/\ 

-T 2 

V 

1.5 
1 

0.5 





(.1 = 1 






= 
— —1 




o 


(i = -2.1 




o 


\i = -3 ^ 




o 






AO 






% 


@ ©oooogg§§ 





1000 



2000 3000 



4000 



5000 



3.5 



-r 2.5 



1.5 



"M=-1 
-(i = 
(i = 1 
(,i = 2.1 
[i = 3 



VrTTl - " " ' ooooooooooooooooo 
Q oo8ooooooooooooo° oouu 

1000 2000 3000 4000 5000 



FIG. 8: MFPT (ti) computed through Eq. (46) as a function of the desorption rate A for several values of the drift coefficient 
for R c = (left) and R c = y/2 > R = 1 (right), in 2D. When \i > 0, the velocity field points towards the origin, while \i < 
means that the velocity field points towards the exterior. For R c < R (resp. R c > R), for a fixed A the search is on average 
faster as \i is more negative (resp. positive). Here d = 2, R = 1, D\ — 1, D2 = 5, e = 0, /c = 00 and a = 0.05 for R c = and 
a = —0.05 for R c — y2 (the truncation size is iV = 200). 
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FIG. 9: The critical ratio Dzc/Di as a function of the target size e in 2D (left) and 3D (right) in the presence of a 1/r velocity 
field with three force intensities: \i — (solid line), \i — 1 (dashed line) and \i — 1.5 (dash-dotted line). The other parameters 
are: R = 1, a = 0.01, R c = and k = 00 (the truncation size is N = 200). 



and V- n {6) = V n {&). These basis functions satisfy 

- A e V n (0) = (tt/0) 2 Pn V n (0) (0 < < 0, n > 0), 
r 2 (V + ^ r ) / n (r) = p\ n \f n {r) (n e Z). 

As previously, we define two scalar products 

(f,g) (f\g)= f(e)g(e)dn d (e), 

JO 

(f,9) (f\g)e= [ f(6)g(6)dfi d (6), 



where d/id(9) is the measure in polar (d = 2) or spherical (d = 3) coordinates for all 9 € [0, 

d m (e) = ^- and dA* 3 (0) = ?^d0- 



(63) 
(64) 



This modified measure is such that the eigenvectors V n {9) are orthonormal: (V n (9)\V n ' (0)} = 5„ 
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Circular sector 



One can easily extend the function g e (0) for a sector of angle . 

g e (6)= 1 -(6-e)(2<J> 

The direct computation yields 

<fle|l> 
(9e\V n ) 

and 

Inn ~ 1 

Inm, 



6). 



(65) 



30 



(<j) — e) cos(7rne/(/>) + 



(j) smfjrne / 4>) 
7r n 



sin(27rne/0) 
27rn 



2m 2 



7r(m 2 



sin(7rne/0) 
cos(7rme/0) cos(7rne/0) 



sin(7rrae/(/>) 



(n > 1), 

(m 7^ n, m, n > 1) 



that generalize formulas from Table [TJ 

In order to complete the formulas for search times, one needs to compute the coefficient 77^ in Eq. (24) and 
the coefficients X n in Eq. (25) that incorporate the radial dependences (e.g., the velocity field v(r) or the partial 
adsorption on the boundary). Since the functions f(r) and fo(r) remain unchanged (see Table \n\ , the coefficient rjd 
is given by previous explicit formulas: Eqs. (51 1 52) with no bias (V = 0) and Eqs. (58, 59) for the velocity field 1/r. 
In turn, the functions f n (r) are modified for the sector because of the prefactor (tt/c/)) 2 in Eq. (64). For instance, if 
there is no bias, / n (0 = r n7r /^, from which 



X n — 



mr /<fi 
kR 



+ (R c /R) 2n7r ^ (x-™^ 



ntr I (f> 
kR 



1 



U7T / (j) 

~k~R~ 



(R c /Ry™/* (1 - ^f) 



that extends Eq. ( |53| ) in 2D. The case of the velocity field 1/r can be studied in a similar way. 
Note that the small e expansion (42) is modified as 



oj 2 T 



2 W 2 (^/7r) 4 x: 



2 W 2 (0/tt) 2 ^ 



X n 



+ 0(e 3 ). (66) 



2. Spherical sector 



One can also compute the MFPT for a spherical sector of angle (j). The angular basis functions V n {6) were given in 
Eq. (62), while the function g e (0) satisfying Eq. (27) with g r e {4>) = is 



... 1 — cos</> f 1 — cos#\ 1+COS0. /l + cos# 
9e{0) = In ( I H In 



1 — cos e 



1 + cos e 



The integration yields 



(1 — COS 0)" , / Sill ( 



In 



2 \ sin e 

One also needs to compute the projections 



1 + cost 6) 2 /l + cose\ cos 6 — cose 

m ( 

1 + cos ( 



(67) 



(68) 



(g e \V n ) e = \/2nTT^ Jd9 sine 5e (0)P„(cos(0 



7T/0)). 



(69) 
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When m = ir/cf) is an integer, cos(m#) can be expressed in powers of cos#, 



[ ^ ] fm-2- i\ 2 m ~ 2 ^- 1 

cos(m#) = 2 m - 1 [cos^] m + m ^ f J 1 : [costf]™" 2 ', (70) 

j=± \ 3 J 3 

(here [m/2] is the integer part of m/2, and we used the convention for binomial coefficients that (™) = 1 for any n) 
so that the computation is reduced to the integrals 

J k = 2k J dOsinO # e (<9)[cos(<9)] fe_1 

e 

,\/r nk ,M Z' 1 _ COSe\ , [cOSe] J ' - [COS (j)V 

= (1 - cos<A)( cos<Ar - l)ln - (1 - cos0) > — 

/, j\/r nfc / /l+COSe\ ,xV^/ , x fc-<i [ COS e ] J _ [ COS 4>] j 

+ (1 + COS0 cos^f - -1 fe In — - l + cos0 V -1 ) k ' [ - ^—^ ^. 

V1+COS0/ J 



3 = 1 



Using this formula, the projections (g e \V n ) e can be easily and rapidly computed. Similarly, one can proceed with the 
computation of the matrix elements I e (n,n'), 

I e {n,n f ) = V^n + lV^n' + 1^ J dOsinO P n (cos (m(9) ) [P n , (cos (m<9)) - P n /(cos(me))] , 

e 

which are reduced to integrals of polynomials. When n/cj) is not integer, the above integrals can be computed 
numerically. 

The radial functions f(r) and /o(r) remain unchanged, while / n (r) are given in Table for the case with no bias. 
The coefficient rjd remains unchanged (cf. Eq. (52)), while the coefficients X n are given by Eq. (60) with 70 = —1/2 



and 7 n = y/n(n + 1)(tt/(/>) 2 + 1/4. The case of the velocity field 1/r can be studied in a similar way. 

3. Multiple targets on the circle 
The MFPT to reach a target of angular extension 2e in a circular sector of half aperture <fi = ir/N t > e (with integer 



m) (see Fig. 10) can actually be rephrased as the unconditional mean search time of N t equally spaced targets of the 
same size 2e on the circle of radius R. Indeed in 2D, due to the reflection principle for random walks, the time spent 
to reach any of the m equally spaced targets on the circle is equal to the time required to reach a single target within 
a wedge with reflecting edges at = ±7r/N t . 



Figure 10 shows the MFPT (ti) in 2D as a function of the number of targets iV t , with a fixed total target length 
^tot = 0.01. This time decreases as l/A/" t 2 , as one can expect from the limiting case A = 0. 

The same procedure in 3D would be to match the time spent to reach any of m equally spaced target caps of size 
2e < lnjNt on a sphere with the time required to reach the target cap G [— e, e] of a cone with reflecting edges at 
— ±7r/7V £ > e (for all (j) G [0, 2ir]). Although not exact because the volume of a sphere cannot be filled by cones, 
this procedure is expected to provide an accurate approximation for the unconditional MFPT as soon as the number 
of targets is sufficiently high. For instance, in the case of 60 equally spaced targets on the sphere, the total excluded 
volume (i.e. the volume between cones) represents less than 1% of the total sphere volume. Knowing that the number 



of membranes or nuclear pores in a cell usually exceeds 100 [12], the results of Sec. IV C should to be relevant for cell 
trafficking studies. 



V. CONCLUSION 



We have developed a general theoretical approach to investigate searching of targets on the boundary of a confining 
medium by surface-mediated diffusion when the phases of bulk and surface diffusion are alternating. This is a 
significant extension of the previous results from [T[ |33] in order to take into account imperfect adsorption, the 
presence of an exterior radial force, multiple regularly spaced targets and general annulus shapes. The coupled PDEs 
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FIG. 10: Left - The search problem of four regularly spaced targets can be represented as a one target search in an angular 
sector 27r/4 with reflecting edges. The shadow green sphere represents the real position of the molecule in the disk while the 
solid green sphere represents its image in the angular sector 2tv/4. Right - MFPT (ti) as a function of the number of targets 
Nt for A = 100, with the total target length e t ot = 0.01. This time decreases as 1/iV* 2 , as one can expect from the limiting case 
A = 0. The other parameters are: R= 1, D± = 1, a = 0.01, R c — 0, k = oo, R c = 0, no bias (V — 0) and D 2 takes three values 
0.5, 1 and 5 (the truncation size is N = 200). 

for the MFPTs t\{6) and ^(r, 6) are reduced to an integral equation for t\(Q) alone whose solution is then found 
in a form of Fourier series. Linear relations for the Fourier coefficients involve an infinite-dimensional matrix whose 
inversion yields an exact but formal solution for the MFPTs. A finite-size truncation of this matrix yields a very 
accurate and rapid numerical solution of the original problem. In addition, we propose a fully explicit approximate 
solution as well as a perturbative one. Although both solutions are derived under the assumption of small targets, 
the approximate solution turned out to be remarkably accurate even for large targets. We illustrate the practical uses 
of the theoretical approach and the properties of the MFPTs by considering in detail several important examples, for 
instance diffusion in a velocity 1/r field. 

The developed approach forms the theoretical ground for a systematic study of surface-mediated processes which 
are relevant for chemical and biochemical reactions in porous catalysts and living cells. From the mathematical point 
of view, the remarkable accuracy of the approximate solution even beyond the expected range of validity remains 
striking and requires further clarifications. 
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Appendix A: Boundary condition for the MFPT 

We check that Eq. Q giving the discontinuity relation of the MFPT between the semi-reflecting surface and 
the bulk can be derived either from a discrete lattice model or from a standard forward equation on conditional 
probabilities [4j [46] . 

1. A discrete lattice approach 

Let us first consider a 2D geometry in which the bulk and surface states are two lattices with radial and angular 
steps Ar and AO. The circular geometry imposes the relation on the radial and angular steps in the bulk at the radius 
r: Ar(r) = rAO. At each time step At, the molecule moves to one of its closest neighboring sites. The value of the 
time step is adjusted according to the position of the molecule: 

At(i,r) = VrAfl/A, 

where i = 1 for the molecule on the adsorbing surface and i = 2 for the molecule in the bulk. This choice maintains in 
the continuous limit a spatially constant value for the diffusion coefficient D2. At r = R, a molecule may either (i) get 
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reflected to r = R — Ar with probability q/2, (ii) get adsorbed onto the surface with probability (1 — q)/2, (iii) move 
along the angular direction, with probability 1/2 (see Fig. 11). Let the random variable r 2 (r,0) (resp. t\{6)) denote 
the first passage time (FPT) for a molecule initially in the bulk at (r, 0) (resp., on the surface at 0). The probability 
for T2(r, 0) to be t = mAt (m G N), is equal to an average of the probabilities of the FPT from neighboring sites to 
be (m - I) At: 

p|r 2 (i?,6>) = mAtJ = |pjr 2 (ii - Ar, 0) = (m - I) At J + !^ipjn(0) = (m - I) At J (Al) 
+ ^\t 2 {R,6 + AO) = (m- l)At| + r 2 (ii, - AO) = (ra- l)Atj. 
The mean FPT in the discrete lattice model is defined as 

oo 

0) = m At p { T 2(r, 0) = mAt}. (A2) 

m=l 



Combining this definition with Eq. (Al) leads to 



t 2 (R, 0) = ±t 2 (R- Ar, 0) + (0) + i (t 2 (R, + A0)+ t 2 (R, - AO)) + At. 



The Taylor expansion of t 2 (R,0) gives 



" a ^,„„,, + T*,„ M = (' - " <"<"> - + (A3) 



Following [H)J 07] , the adsorption coefficient is denned as 

, _ l-q 1 



9 i?A6>' 



(A4) 



In the continuous limit, when all Ar, AO, At tend to with Ar/AO = r and D 2 = Ar 2 /(2At) constant, Eq. (A3) 
turns out to be expressed in terms of k only: 

^ =k{t 1 (6)-t 2 (R,6))+o(l). (A5) 

ar | r =(i?,6>) 

Indeed At/ Ar = o(l) and AO 2 / Ar = o(l) in this limit. For a perfectly adsorbing boundary we have q = 0, k — >• oo, 
which is indeed compatible with the continuity relation t 2 (R,0) = used in Ref [1 . For a perfectly reflecting 

boundary g = 1, k — » and d r t 2 (r, 0) = at r = [4]. This is indeed the condition that we imposed on the boundary 
r = i? c in Sec. 11111 




FIG. 11: A discrete lattice model. Regular AO slicing imposes At and Ar to vary over the domain S in order to maintain an 
isotropic diffusion constant D 2 . Only one quadrant is presented. 
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2. Equivalence with a forward boundary condition 



We now check that the MFPT condition of Eq. Q is also compatible with the following boundary condition on 
the conditional probability [4] 



dp((r,0),t\x,t') 

dr \ r =R 



kp((r,6)\x,t% 



\r=R, 



(A6) 



where p(x, t\x f , t r ) is the probability for a molecule to be at x at time t provided that the molecule was at x' at an 
earlier time t' < t. We denote the spatial coordinate x = (r, 6) if the molecule is in the bulk and x = 6 if it is adsorbed 
on the surface. 

We follow the standard method presented in [46 j. The stochastic process under study is Markovian hence the 
conditional probabilities satisfy the Chapman-Kolmogorov equation, with t > s > t' 

p(x, t\x\ t') = [ du(y) p(x, t\y, s)p(y, s\x f t') + / R^d^O) p(x, t\0, s)p(0, s\x\ t'), 
Js Jo 

where S = (R, R c ) x [0, it] and the measure dv is 

di/ 2 (r,0) = 2 r drdO, dv 2 (6) = 2 d0, di/ 3 (r,0) = 27rrsin(9 drdO, dv 3 (0) = 2tt sin OdO. 

Taking the derivative with respect to the intermediate time s leads to the expression 



0=j- s p(x,t\x',t') = Jdu(y) dp{x fy> s) p{y,s\x>,t>) 



7 / s\\ dp(x.t\6.s) , , ,s 
dv{9) PK g J J p(0,s\x',t') 



dv(y) p{x,t\y,s) 
dv(0) p(x,t\6,s) 



dp(y,s\x',t') 
ds 

dp(e,s\x',t') 
ds 



The backward Chapman-Kolmogorov equations read 
dp{x,t\6,s) _ D 1 



ds 

dp(x,t\(r,9),s) 
ds 



A e p(x, t\0, s) + A<^ p(x, t|0, s) - p(x, t\R - a, 0, s) 



R 2 



= -D 2 A M) p(x,t\(r,9),s)-v(r) Vp{x,t\(r,0),s). 



(A7) 

(A8) 
(A9) 



The forward Chapman-Kolmogorov equations are 
dp(0,s\x',t') 



ds 

dp(y,s\x',t') 
ds 



^A 6P (0, s\x',t') - A p(6, s\x>, «) - D 2 ^ *>' ^' 



dr 



\r=R 



+ v(R) p((R,0),s\x',t'), 
D 2 A y p(y, s\x', t') - V (v(r)p(y, s\x', t')) 

d-l 



+ A 



R 



R — a 



5 d (r-(R-a,6)) p(6,s\x' ,t'). 



(A10) 



(All) 



The terms in Eqs. (A10) and (All) are justified as follows: (i) —A p(0, s\x' , t') corresponds to a constant rate of 
desorption from the surface to the bulk; (ii) —D2d r p((r,Q),s\x',t') is the flux into the surface due to diffusion; (hi) 
v(R) p{{R,6), s\x r ,t f ) is the flux into the surface due to the drift (by convention, v(R) > for a velocity drift field 
oriented towards to the exterior); (iv) A [R/(R — a)] _1 6 d (r — (R — a, 6)) p(0,s\x' ,t') corresponds to the flux into 
the bulk due to the desorption from the surface and the ejection at a distance a (5 being the Dirac delta function). 

For convenience, we will use the shorthand n otat io ns p( y, s\x' , t r ) = p( y), p(x, t\y, s) = p(y) and p(6, s\x' , t r ) =p(0). 
Substituting the Chapman-Kolmogorov Eqs. (A8 - |A11| ) into Eq. ( |A7[ ) leads to the following equation 







dv(y) [-D 2 A yeS p(y)]p(y) 
dv{y) p(y) 

[ Rd ~ ldm & 



D 2 A ye s p(y) + A 



R 



d-l 



R-a 

A e p(0) + \{p(0,s)-p(R-a,0)} 



S d (r-(R-a,0)) p{0) 



p(0,s) 



f 

Jo 



R d ~ l dv(0) p(6) 



—^Aep(O) - A p{6) - D^—r 

R d or \r=R 
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One notices that 



/ du(y) 
Js 



R 



R — a 



d-l 



S d (r -(R-a, 0))p(r, 0) p(0) + / R d ~ x dv(0)p(R - a, 0)p(0) = 0, 

Jo 



and that the four terms proportional to A vanish. Two terms with angular Laplacians also cancel each other due to 
the hermit icity of the angular diffusion operator: 



f dv{0) A 6 p(6)p(y) - f dv{9) A p(0)p(y) = 0. 
Jo Jo 



The divergence theorem yields the integral over the frontier dS of the annulus S: 



= D 2 I dv{6) 

IdS 



dp(r,0) dp(r,0) dp(r,0) 
p(r,0) p(r,0)- 



dr \r=R 



This equality can be satisfied only if: 
dp((r,9),s) 

dr \r=R 



p((R,6),s) 



dr \r=R 

dp((r,0),s) 

dr \r=R 



dr \r=R 



p((r,e),s)-p(8,s) 



p(0,s) 



(A12) 



Inserting the forward boundary condition (A6) into Eq. (A12) gives the boundary condition on the backward proba- 
bility distribution 



dp(x,t\(r,6),s) 

dr \ r =R 



= k 



p(x,t\0,s) -p(x,t\(r,0),s) 



\r=R 



Integrating over the space and time variables x and t, we obtain the boundary condition for the MFPT: 



dt, 



=k{t 1 (0)-t 2 (R,0)} (0<0<7r), 

Or \r=(R,G) 



which identifies with Eq. Q. 



Appendix B: Interpretation of rjd/D 2 as a mean first passage time 

We consider the probability density H(0\6) for a molecule initially at the bulk point (R — a,Q) to first reach the 
surface r = R at the angle 0. The mean duration of this Brownian path is denoted t c (0\0). 

The MFPT t 2 (R — a,Q) to reach the target can be expressed as the averaged sum of the MFPT to reach a point 
(i?, 0) on the surface and the MFPT to reach the target from this point of the surface, the probability density for the 
first hitting point (R,0) being the harmonic measure H(0\0) : 



t 2 (R - a, 0) = J (t c (0\0) + U(0\0)dfji d 



(9). 



In 2D and in the general case considered in Sec. Ill, the probability density 11(0 10) is 

oo 

II(0|0) = 1 + 2 ^(X„ + 1) cos(n(£ - 9)), 



(Bl) 



(B2) 



where X n is given by Eq. (25). Substitution of this expression in Eq. (Bl) leads to 



t 2 (R-a,0) = (t 1 ) + - t c (6\6)IL(6\6)d0 + - + 1) / cos(n((9 - O^t^S. 

n Jo * n=1 Jo 



Identification with Eq. ( 36 ) gives 



D 2 



- I t c (e\0)n(0\e)de = - [ t c (9\o)ii(9\o)$, 

n Jo n Jo 



(B3) 



(B4) 
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which identifies '-D2 as the MFPT to the circle of radius i?. In particular, it can be shown that in the 2D case of 
Sec. II V A II 

U(0\0)t c (0\0) = ^ (1 - (r/i?) 2 ) ^1 + fj cos(n(0 - 0))^ . (B5) 



One can verify that the substitution of this expression into Eq. (|B4j) leads to the well known result of Eq. (|47|), 
m = R 2 (l-x 2 )/A. 



The argument leading to Eq. ( B4 ) can be extended to the 3D case with the following expression for the probability 
density: 

00 

II(0|0) = 1 + ^(2n+ 1)(X„ + l)P„(cos0)P n (cos0). (B6) 



Appendix C: Matrix elements I e (n,m) in 3D 

The matrix elements 7 e (n,m) in 3D were computed in [1 J. An explicit formula for non-diagonal elements (m ^ n) 
is given in Table [TJ In turn, the diagonal elements I e (n,n) can be expressed as 

I e (n,n) = -P n (u) — + — — — — , u = cose, 

n + 1 2n + 1 

through the function F n (u), for which the explicit representation was derived in [1 
F n (u) = u[Pl{u) + 2P*_ 1 (u) + ... + 2Pl{u) + P (u)} 



- 2P n (u)P n _ 1 (u) - 2P„_ 1 (u)P„_ 2 (u) - ... - 2P 1 (u)P (u) + u 

= £ [2(u - l)Pi(u) + [P k (u) - P fc _!( W )] 2 ] - (u - 1)P» + (u - 1)P 2 ( W ) + u. 



(CI) 



k=l 

One can also check that this function satisfies the recurrence relations 

F n (u) = F n -i(u) + u[P 2 (u) + P n 2 _!W] - 2P n (u)P n - 1 (u), F (u) = u. (C2) 
that simplifies its numerical computation. Note that F n (=bl) = F n _i(=bl) = ... = ±1. 

Appendix D: Case of a 1/r 2 velocity field 

We now examine the 3D case of a radial 1/r 2 velocity field v(r), which is characterized by the dimensionless 
parameter fi: 



The function / is expressed as 



fiD 2 R _> , v 

v[r) = v. — r. (Dl) 



/» = - r l - r *± + ^ e -^Ei(^/r), (D2) 

o 



where Ei(z) is the exponential integral: 

z 

f e x 

Ei(z) = / — dx. 

-co 

The function fo is 
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(this particular choice of the additive and multiplicative constants ensures that R/r is retrieved in the limit /i — )> 0). 
Radial functions f n (r) are found as products of powers and confluent hypergeometric functions 1F1 of r: 



fn(r) 
f-n(r) 



iF^-n, -2n, -fiR/r) 



(n > 0), 



iFi(n+l, 2n + 2, -jiR/r) (n > 0). 



(D4) 
(D5) 



For n = 0, this expression is reduced to e ^ R l r . In the limit \i — » 0, the above functions reduce to r n and r n 1 from 
the earlier case /i = 0. On the one hand we have 

r / n (r) = w^-SFit-n+l, -2n, -/xii/r) (n > 0), 

d r f-n(r) = -(n + l)r- n - 2 iF!(n + 2, 2n + 2, -/xii/r) (n > 0), 



from which 



On the other hand we have 



d r fn(r) 



d r f- n (r) n+1 



r 2n+i + -2n, -fiR/r) 

iFi(n + 2, 2n + 2, -fiR/r)' 



from which 



d r fo(r) 
dj(r) _ 



^V^Ei W) - ^ - - 2r 



i? 



-nR/r 



(pR) 2 Ei(nR/r) - e^ R/r [fxRr + r 2 - 2r 3 /(jiR)] 



d r fo(r) 

This last expression is needed to compute the quantities % and X n . 



(D6) 



(D7) 
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